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Abs tr act 

The present paper studies the numerical simulation of flows 
with shock/boundary-layer upstream interaction, under con- 
ditions of symmetry in geometry, boundary conditions, and 
grid. For this purpose, a series of two- and three-dimensional 
numerical test-cases were carried out. The tests showed that 
standard numerical schemes, which appear to be symmetry- 
preserving under most flow configurations, produce non- 
symmetric perturbations when large separated regions are 
present. These perturbations are amplified when the core 
flow is under compression. If the flow-blockage due to sepa- 
ration is sufficiently large, the symmetry of the flow may col- 
lapse altogether. Experimental evidence of this numerical 
behavior is also considered. 
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Nomenclature 

H test-section inlet height [ml 

k turbulence kinetic energy [m 2 /s 2 ] 

M Mach number 

p pressure [Pa] 

T temperature [°K] 

V velocity [m/s] 

x axial (streamwise) coordinate [m] 

y vertical coordinate [m] 

y + dimensionless vertical turbulence coordinate 

z lateral coordinate [m] 

|i viscosity [N s/m 2 ] 

p mass density [kg/m 3 ] 

co specific dissipation rate [ 1/s] 

Subscripts: 

0 total conditions 

back back (exit) conditions 

exit exit conditions 

inlet inlet conditions 

j injectant (jet) conditions 

s static conditions 


L Introduction 

In a previous paper 1 the author and co-workers 
reported an unexpected numerical behavior in the simula- 
tion of dual-mode scramjet combustors. Under conditions 
of perfect symmetry in geometry, grid and boundary con- 
ditions, computational fluid dynamic (CFD) calculations 
showed an asymmetric response when a shock-train/sepa- 
rated-flow system was present in the flowfield. Ollier 
researchers using a completely different computational 
methodology reported similar results 2 . 

The problem under consideration is shown in figure 1 ? 
which depicts the transverse plane of a flow in a rectangu- 
lar cross-section duct. A supersonic incoming flow is sub- 
jected to an adverse pressure-gradient due to a back-pres- 
sure Pexit» the source of this pressure could be, for 
example, a combustion process or a diffuser. In general, 
for any p exit between inlet and normal-shock values the 
flow configuration will be a confined, compressed core- 
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flow surrounded by separated regions 2 3 . The separating 
boundary-layer creates a series of shocks, or shock train, that 
further compresses the core flow. This shock-train will be 
mainly oblique or normal, depending on whether the exit 
Mach is supersonic or subsonic. Since the sbock-train/sepa- 
rated flow interaction propagates upstream of the source of 
pressure, the whole phenomenon is referred to as ‘upstream 
interaction’. The symmetry and stability of flows with 
upstream interactions is the focus of the present paper. 

There exists experimental evidence to suggest that these 
flows may not remain symmetric. Carroll 4 and Carroll and 
Dutton 5 reported a series of experiments on supersonic rect- 
angular ducts subjected to adverse pressure-gradients. They 
were able to create an upstream interaction similar to the one 
depicted in figure 1. At low Mach numbers (M - 1.60) the 
flow appeared to be mostly symmetric. At high Mach num- 
bers (M ~ 2.50), however, the schlieren visualization and 
pressure traces showed a large degree of asymmetry, similar 
to the one reported in reference 1 

Asymmetrical behavior under symmetry conditions 
have been previously reported for other types of flow. Levy 
et al 6 did so for flows around slender bodies at high angles of 
attack (AOA). Experimentally, these flows show a convective 
instability mechanism, i.e., minute imperfections can cause 
large flow asymmetries which vanish when the disturbance is 
removed. The authors contended that “numerical simulation 
of such flows should give the basic symmetric solution as 
long as algorithms, boundary and initial conditions, and 
grids are perfectly symmetric”. They showed that the left- 
hand-side of some implicit time-integration algorithms could 
create asymmetries in the steady-state solution. In particular, 
the diagonal approximate-factorization (DAF) algorithm of 
Pulliam and Chaussee 7 was shown not to be symmetry-pre- 
serving. For low AOA the effects were small, the flowfield 
appearing largely symmetric. At high AOA, however, the 
flowfield was manifestly asymmetric. Furthermore, the 
asymmetry did not correlate with the experimentally- 
observed pattern. 

There may be some analogy between external flows 
around high-AOA slender bodies, and internal flows with 
upstream interaction. The duct flows at low Mach numbers 
would resemble the low-AOA bodies, with the flow largely 
symmetric; the asymmetry of high-Mach duct-flows would 
be similar to that of high-AOA bodies. Unlike the slender- 
body experiment, the asymmetry in the duct flows could not 
be traced to any known perturbation. This would indicate a 
more fundamental instability than the convective type of the 
slender-body flows. 

The present paper will attempt to determine the nature 


of the numerical asymmetry observed in flows with 
upstream interaction. First, a series of two-dimensional 
numerical tests will be carried out to determine the possi- 
bility of obtaining a basic symmetric solution for flows 
with upstream interaction. Next, diree-dimensional effects 
will be studied in a rectangular-section duct. Finally, the 
experiment of Carroll and Dutton will be modeled to com- 
pare numerical asymmetries with experimental ones. 

2. Computatio nal Approach 

All calculations presented in this paper were done 
with NASA Langley’s VULCAN code. This code has been 
extensively described before 8 . Only the specific implemen- 
tation used in this paper will be detailed here. 

The inviscid fluxes were modeled with Edwards low- 
dissipation flux-split scheme 9 . Fluxes were evaluated 
either with a first-order approximation, or with a third- 
order MUSCL scheme and Van Leer’s flux limiter. 

The full viscous fluxes were modeled with second- 
order central differences. In all cases, viscous flows were 
assumed fully-turbulent. Wilcox’ 1998 k-co turbulence 
model was used, with Wilcox’ wall-matching functions at 
solid walls 10 . 

In all cases, the flow was assumed to be spatially- 
elliptic. The gas was modeled as calorically-perfect air. 
Time-integration was done with either a two-step explicit 
Runge-Kutta (RK) scheme, or the implicit diagonal 
approximate-factorization (DAF). 

All calculations were run in double-precision on 64- 
bil machines (Cray c90. Origin 2000, SGI R 10000) 


3. Two-Dimensional Nu merical Tests 

A series of two-dimensional numerical test cases were 
devised to establish the existence of a numerical scheme 
that could produce a basic symmetric solution for flows 
with upstream interaction. Each case would try to isolate 
or add elements present in these flows, and determine their 
effect on the symmetry of the solution. The tests were 
based on the conditions for the Carroll and Dutton 
experiment 5 . The inlet conditions (assumed uniform for all 
cases) were M = 2.50, pp = 3 10 KPa and Tq = 295 °K The 
inlet height was H = 0.0381 m (used to non-dimensional- 
ize lengths in all plots for this section). 

In all cases, inviscid fluxes were modeled as first- 
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order, in order to find the simplest algorithm that would be 
symmetry-preserving. Unless otherwise noted, all wall con- 
ditions were no-slip adiabatic and supersonic extrapolation 
was imposed at the exit. 

34 1 Co nv e rgent Vw t 

Figure 2 shows a convergent duct, which is reminiscent 
of the core flow in figure 1 . The main purpose of the test is to 
produce a shock train similar to the one observed in 
upstream interaction but without any separated flow being 
present. This would allow a preliminary determination on the 
symmetry-preserving properties of a given algorithm. The 
grid is symmetric with respect to y = 0. Its dimensions are 
177 x 81 (stream wise x height). All results to be shown for 
this case were considered converged after the L 2 -norm of the 
residual decreased four orders of magnitude. 

The Mach contours obtained using RK integration and 
assuming inviscid flow are presented in figure 3-a); a series 
of shocks can be seen, similar to the shock train in the core- 
flow of figure 1. The flow appears symmetric and this is con- 
firmed by the streamwise variation of the vertical velocity 
component at the symmetry plane, V y (y = 0) (figure 3-b): 
the distribution is mostly noise, without any apparent corre- 
lation with the flowfield. Adding the viscous fluxes (figure 4) 
produces a qualitatively similar flowfield; the resulting 
scheme also appears to be symmetry-preserving. 

Integrating the inviscid flow with the DAF scheme 
shows a Mach-number pattern virtually identical to that of 
the RK method (figure 5-a), but now V y (y = 0) is no longer 
noise (figure 5-b). Its distribution is several orders of magni- 
tude bigger than that of the RK scheme, and it correlates 
with the flowfield; in particular, the two peaks correspond to 
the shocks intersecting the symmetry line. This behavior of 
the DAF scheme agrees with the findings of Levy et al 6 . 

This simple test seems sensitive enough to capture the 
symmetry-preserving behavior of an algorithm, and suggests 
that at least first-order space-integration with the RK time- 
integration is symmetry preserving. It would also indicate 
that, by itself, a shock-train is not a fundamental source of 
asymmetry. 

3,2, Sudden-Expansion Duct 

Figure 6 shows the layout of a duct with a sudden 
expansion or backward-facing step. This domain wifi be 
used to generate a shock-train/boundary-layer interaction. 
This can be accomplished two ways: by injecting air at a 
pressure p in j and very low speed (as shown in the figure), or 
by applying a back-pressure Pbzck at the exit. Either proce- 


dure can create a shock/separated-flow interaction, 
whether upstream or not. Injection makes it easier to con- 
trol the interaction at low or moderate compression in the 
core flow; back-pressure is better suited for high compres- 
sion. As will be seen, the degree of core -flow compression 
will have an effect on the symmetry of the flow. Unless 
otherwise noted the grid dimensions are 153 x 161 (inlet 
or isolator) and 201 x 201 (expansion). 

Figure 7 presents the results for p inj = 40 KPa, using 
the RK scheme. The isolator has been shortened to reduce 
computational effort and its grid has dimensions 53 x 161 . 
These calculations were converged six orders of magni- 
tude. The Mach contour-plot (figure 7-a) shows an oblique 
shock-train interacting with the separated flow produced 
by the injection. The core-flow exit-pressure is slightly 
larger than the inlet pressure (about 20 KPa). The flowfield 
appears symmetric, but the centerline normal-velocity dis- 
tribution (figure 7-a) is greater than machine-zero and cor- 
related to the flow. The results obtained using the DAF 
scheme are virtually identical (figure 8), except for the fact 
that the centerline perturbations are an order of magnitude 
bigger. 

Increasing the injectant pressure to p inj = 95 KPa pro- 
duces substantial upstream interaction, as can be seen in 
figure 9 (after converging six orders). The core-flow exit 
pressure turns out to be twice as large as the inlet pressure, 
indicating a higher core-flow compression with respect to 
the previous case. Apparently as a result, the velocity per- 
turbation generated with the RK scheme has increased by 
two orders of magnitude. When using the DAF scheme 
(figure 10), the flow symmetry collapses altogether. 

Figure 11 shows the results (converged four orders) 
when the upstream interaction is brought about by apply- 
ing a back-pressure at the exit; in this case, p^^ = 80 KPa. 
or about four times the inlet pressure (and approximately 
twice the core-flow compression of the previous case). As 
can be seen, the flow symmetry collapses regardless of die 
numerical scheme (RK or DAF) used. For comparison, fig- 
ure 12 shows what the flowfield would look like if only 
half of the duct height had been solved, with a symmetry 
assumption imposed at the centerline. Not only the flow is 
symmetric, but the extent of the upstream interaction into 
the isolator is about half as much as in die previous case. 

3,3. Analysis 

From the previous series of numerical tests, some pre- 
liminary conclusions may be drawn: 

• The presence of large separated-flow regions interacting 
with shock-trains appears to be the immediate source of 


3 

American Institute of Aeronautics and Astronautics 



asymmetric perturbations greater than machine-zero. 

• These perturbations are present even with numerical 
schemes that appear to be symmetry-preserving when 
applied to more conventional flows. 

• The asymmetries are either introduced numerically (by the 
scheme or the modelization itself), or are inherent to flows 
with large shock/separation interaction. 

• If the latter is true, a basic symmetric solution (i.e., sym- 
metric within machine zero) for these flows may not be pos- 
sible with conventional computational techniques 

• Regardless of their origin, asymmetric perturbations seem 
to increase with the core-flow compression. 

• Flows that are that are under relatively high compression 
appear to be unstable, i.e., any perturbation will be amplified 
until the flow-symmetry collapses regardless of the scheme 
used. 

• The extent of the upstream interaction appears to be larger 
when the flow symmetry is allowed to collapse, than when 
the symmetry is imposed as boundary condition. 

4. Three- Dimensional Numerical Test 

To determine the numerical stability of three-dimen- 
sional flows, a duct with rectangular cross-section was 
tested; it is based on the isolator section of the experimental 
scramjet modeled in reference 1. Figure 13 shows the layout 
of the duct, with the two symmetry planes also drawn. The 
duct height was H = 0.032 m (used to non-dimensionalize all 
plots in this section). Inlet nominal conditions were M = 
2.50, p 0 = 1.0 MPa and T 0 = 2000 °K. To promote the exist- 
ence of upstream interaction, an incoming boundary-layer 
profile was imposed at the inlet. A step-wise uniform back- 
pressure distribution was applied at the exit: a value of 150 
KPa was used a small distance from the sidewalls, and 250 
KPa elsewhere in the section. It was hoped this distribution 
would recreate the expected flowfield of a dual-mode scram- 
jet. For qualitative purposes a coarse grid of dimensions 101 
x 25 x 125 (streamwise x height x width) was used. The 
numerical scheme consisted of third-order fluxes with DAF 
time integration. Since the flow is highly compressed (exit 
pressure was four times the inlet for the most part), a sym- 
metry-preserving algorithm was not deemed necessary. 

The overall flowfield configuration can be seen in figure 
14. A central separated-flow region, similar to the two- 
dimensional case, is clearly visible. Separation adjacent to 
the sidewalls is also apparent. The Mach contour-plot in the 
vertical symmetry-plane (figure 15-a) shows that flow sym- 
metry has collapsed, just like in the two-dimensional case. In 
the horizontal symmetry-plane (figure 15-b), symmetry has 


also been lost, but not to the same extent as in die vertical 
plane. This may be due to the fact that flow blockage from 
separation is smaller in the horizontal plane titan in die 
vertical. If this is the case, it may suggest diat die stability 
of the flow-symmetry depends, for a given amount of com- 
pression, on the ‘slenderness’ of die core flow widt respect 
to the geometric duct: the narrower the core-flow relative 
to the relevant duct dimension, the more likely the flow 
will depart from symmetry or its symmetry collapse alto- 
gether. 

5. Experimental Test Case 

The results from the previous secuons showed dial 
asymmetry is associated with flows containing upstream 
interaction, and with highly-compressed, slender core- 
flows. This asymmetry will now be compared widt die one 
observed experimentally in the series of experiments car- 
ried out by Carroll and Dutton 4 5 . Their full description may 
be found in the references cited previously; only a brief 
oudine will be presented here, followed by die numerical 
results. 

5.1. Descrip tion of the Experiment 

The schematic of the supersonic wind-tunnel facility 
is shown in figure 16. The wind tunnel can be divided in 
four sections, in downstream order: a constant-area stagna- 
tion chamber; a convergent-divergent nozzle; a test-section 
with divergent walls; and an exit diffuser. The entire facil- 
ity is symmetric, and the cross-section, rectangular widi 
constant width. For the nominal test M inle , = 2.50, the inlet 
height is H = 0.0381 (which again was used for non- 
dimensionalization of the results); the section aspect-ratio 
is 2. The top and bottom walls have a divergence angle of 
0.25 deg. The operating conditions are po = 310 KPa and 
T 0 = 295 °K. 

Figure 17 presents experimental oil and schlieren 
visualizations, for a diffuser geometry such that the test- 
section exit back-pressure is p ex „/po = 0.374. The oil visu- 
alizations correspond to the top and bottom surfaces, while 
the schlieren provides a side view. The asymmetry of the 
flow is evident. Following the terminology of the research- 
ers, the flow is said to be “attached” to the lower wall and 
“separated” from the top wall (even though separation is 
apparent in both walls). Note also the three-dimensionality 
of the flow, with comer recirculation possibly due to the 
presence of secondary flows. There appears to be a slight 
asymmetry in the horizontal plane, especially on the top 
surface. 
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$t2t Num erical Res u lts 

The test was approximated with a two-dimensional 
assumption. The stagnation and nozzle sections were solved 
first, and the exit conditions from the nozzle were used as 
inlet conditions for the test section. Each domain was solved 
elliptically using third-order inviscid fluxes and the DAF 
scheme. Since it is expected for this flow to be unstable, the 
DAF scheme was chosen to allow high CFL numbers. Con- 
vergence was assumed when the wall-pressure traces 
remained essentially unchanged. 

The grid used for the stagnation and nozzle domain is 
shown in figure 18, with an expanded view of the nozzle. 
Only half of the domain was solved, with a grid of dimen- 
sions 435 x 101. The inlet conditions were subsonic, and 
extrapolation was used at the exit. A symmetry condition 
was imposed at the top, and no-slip adiabatic was assumed at 
the wall. For this domain no wall-functions were used, the 
boundary-layer being integrated to the wall; the resulting y + 
had a maximum value of about 1 at the throat. Calculations 
were converged four orders of magnitude. The resulting 
Mach contours can be seen in figure 19. 

The computational domain for the test section extended 
to its entire length, l test /H = 20. The exit diffuser was repre- 
sented by a mixed subsonic/supersonic exit boundary condi- 
tion, with the experimental pressure applied there. The 
baseline grid had dimensions 801 x 161. The grid was clus- 
tered in the streamwise direction between x/H values of 1 
and 7, region where the leading-edge of the shock-train was 

expected to be found. The maximum y + turned out to be 
approximately 40, which required the use of wall-functions. 

Figure 20 shows the Mach contours for the initial half of 
the domain (note the vertical scale is twice the horizontal, for 
purposes of clarity). The asymmetry has taken place, and 
several pairs of oblique shocks are visible. The flow is 
attached to the bottom and separated from the top as in the 
experiment, but this is entirely fortuitous; whether the flow 
separates on one surface or the other is entirely random. 

The wall-pressure plots for the separated wall are plot- 
ted in figure 21. The location of the shock leading-edge is 
close to the experiment. The slope of the pressure rise, how- 
ever, is underpredicted. Several factors may account for the 
disagreement. The most important may be the three-dimen- 
sional effects pointed before; comer effects may constraint 
the core flow, increasing the pressure gradient with respect to 
two-dimensional flows. Replacing the diffuser with a uni- 
fonn back-pressure could also have an effect. 

In the same figure, results are shown for a half-height 


domain with a symmetry boundary-condition at die center- 
line. As before, the extent of die upstream interaction is 
considerably smaller and gives a much-poorer agreement 
with the data. 

The effects of grid resoludon can be seen in figure 22. 
The medium grid corresponds to die baseline grid 
described before. The dimensions of die coarse and line 
grids differ from the medium by a factor of two in each 
coordinate direction. The difference between the medium 
and fine grid appears small enough (and much smaller dian 
between coarse and medium) to suggest grid convergence. 


6. Conclusions 

A series of numerical test cases was conducted to 
study the behavior of flows widi upstream interaction. In 
particular, the possibility of obtaining a basic symmetric 
solution under conditions of symmetry in geometry, 
boundary conditions and grid was examined. Widi die 
usual numerical schemes, it was found impossible to 
obtain a basic symmetric soludon. The results showed dial 
the presence of large separated- flow regions interacting 
with shock-trains create non-symmetric perturbations that 
are greater than machine-zero. These pcrturbadons, when 
coupled with the strong compression characterisdc of 
flows with upstream interaction, may bring about die loss 
of the symmetry in the flowfield. The extent of departure 
from symmetry appears to be a funcdon of die flow-block- 
age produced by the separadon itself. The larger die asym- 
metry, the longer the extent of die upstream interaction 
compared to a symmetry-imposed calculadon. 

If the above is true, it may have an impact in some 
practical applications. Operational configurations, such as 
actual scramjet engines, are characterized by large aspect- 
ratio sections. Therefore, die modeling of half of the sec- 
tion with a symmetry boundary condition at the vertical 
symmetry plane is probably adequate. For small aspect- 
ratio sections (used in most experimental configurations) 
the use of any symmetry assumption is more problemati- 
cal. The critical aspect ratio has yet to be determined. 

An attempt was made to simulate the only experiment 
known to the author that displayed non-symmetric behav- 
ior in a flow with upstream interaction. The results were 
able to capture the qualitative behavior of die experiment. 
Quantitative comparison probably suffered from the use of 
a two-dimensional assumption to model a low aspect-ratio 
three-dimensional flow. A full three-dimensional simula- 
tion may be needed to determine die correlation, if any, 
between physical and numerical asymmetry. From die pre- 
vious conclusions, the full cross-section may have to be 
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modeled. 

It would also be instructive to study the stability of flows 
with upstream interaction within annular ducts. Experimen- 
tal evidence" 12 suggests that geometrically-axisymmetric 
flows remain fairly symmetric in the presence of upstream 
interaction. Whether their three-dimensional, full-section 
simulation also displays any symmetry remains to be deter- 
mined. 
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b) Vertical velocity at y = 0. 



Shock train 



Figure 1: Schematic of a rectangular duct with 
shock-train/boundary-layer interaction. 


Figure 3: Convergent duct: Inviscid flow with 
RK integration. 


a)Mach Contours 
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Figure 2: Convergent duct: layout. 


a)Mach Contours 



b) Vertical velocity at y = 0. 



Figure 4: Convergent duct: Viscous flow with 
RK integration. 
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a)Mach Contours 


a)Mach Contours 



b) Vertical velocity at y = 0. 



Figure 5: Convergent duct: Inviscid flow with 
DAF integration. 


b) Vertical velocity at y = 0. 



Figure 7: Sudden-expansion duct: 
Pinj = 40 KPa - RK integration. 


a)Mach Contours 




Figure 6: Sudden-expansion duct: layout. 
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v, (y = 0) 


b) Vertical velocity at y = 0. 



Figure 8: Sudden-expansion duct: 
Pinj = 40 KPa - DAF integration. 
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Figure 10: Sudden-expansion duct 
Pinj = 95 KPa - DAF integration. 


a)Mach Contours 


a) RK 
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b) Vertical velocity at y = 0. 



b) DAF 


Mach 
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Figure 9: Sudden-expansion duct: 
Pinj = 95 KPa - RK integration. 


Figure 11: Sudden-expansion duct: 
Pback = 80 KPa. 
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a) Vertical symmetry.-plane 
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Figure 12: Sudden-expansion duct: 
Pback = 80 KPa - Symmetry assumption. 
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b) Horizontal symmetry-plane 


Mach 




Figure 15: 3-D isolator: Symmetry planes. 


Figure 13: 3-D isolator: layout. 


Mach 



Figure 14: 3-D isolator: Overall flowfield. 



Figure 16: Experimental test-case: 
Layout of the facility (courtesy B.F. Carroll). 
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Figure 17: Experimental test-case: Figure 20: Experimental test-case: 

Oil and schlieren visualizations Test-section Mach contours, 

(courtesy B.F. Carroll). 



Figure 18: Experimental test-case: 
Stagnation and nozzle section grid. 



Figure 21: Experimental test-case: 
Test-section wall pressures - Baseline grid. 


Mach 



Figure 19: Experimental test-case: 
Stagnation and nozzle section Mach contours. 



Figure 22: Experimental test-case: 
Test-section wall pressures - Grid resolution. 
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